      PROGRAM PGDEM4
C-----------------------------------------------------------------------
C Test program for PGPLOT: test of imaging routine PGGRAY (with 
C PGCONT). This program serves to both demonstrate and test PGGRAY.
C It computes an (arbitrary) function on a 2D array, and uses both
C PGGRAY and PGCONT to display it. An irregular transformation (TR 
C matrix) is used, to test (a) that the routine works when the array
C pixels are not aligned with the device pixels, and (b) that the
C image is clipped correctly at the edge of the viewport. The program
C also draws the bounding quadrilateral of the contour map. The contours
C should end on this quadrilateral, but note that the grayscale image 
C extends one half pixel beyond this quadrilateral (if the subarray to 
C be displayed contains N pixels in one dimension, the width of the
C image is N units, while the width of the contour map is N-1 pixels).
C-----------------------------------------------------------------------
      INTEGER PGBEG
      INTEGER   MXI, MXJ
      PARAMETER (MXI=40, MXJ=40)
      INTEGER I,J
      REAL F(MXI,MXJ)
      REAL ANGLE,FMIN,FMAX,ALEV,RADIUS,TR(6)
C
C Open device for graphics.
C
      IF (PGBEG(0,'?',1,1) .NE. 1) STOP
C
C Compute a suitable function.
C
      FMIN = F(1,1)
      FMAX = F(1,1)
      DO 20 I=1,MXI
          DO 10 J=1,MXJ
              F(I,J) = COS(0.6*SQRT(I*80./MXI)-16.0*J/(3.*MXJ))*
     :                 COS(16.0*I/(3.*MXI))+(I/FLOAT(MXI)-J/FLOAT(MXJ))
              FMIN = MIN(F(I,J),FMIN)
              FMAX = MAX(F(I,J),FMAX)
   10     CONTINUE
   20 CONTINUE
C
C Compute the coordinate transformation matrix.
C ANGLE is the image rotation angle in radians.
C
      ANGLE=-120./57.29578
      RADIUS= 40.*SQRT(2.)
      TR(1) =-(40./MXI)*COS(ANGLE)+(40./MXJ)*SIN(ANGLE)-
     :         RADIUS*COS(ANGLE+45/57.29578)
      TR(2) = 80.*COS(ANGLE)/MXI
      TR(3) =-80.*SIN(ANGLE)/MXJ
      TR(4) =-(40./MXI)*SIN(ANGLE)-(40./MXJ)*COS(ANGLE)-
     :         RADIUS*SIN(ANGLE+45/57.29578)
      TR(5) = 80.*SIN(ANGLE)/MXI
      TR(6) = 80.*COS(ANGLE)/MXJ
C-----------------------------------------------------------------------
C PGGRAY
C-----------------------------------------------------------------------
C Clear the screen. Set up window and viewport.
C
      CALL PGPAGE
      CALL SETVP
      CALL PGWNAD(-40., 40.,-40., 40.)
      CALL PGSCI(1)
C
C Draw the map with PGGRAY.  
C
      CALL PGGRAY(F,MXI,MXJ,1,MXI,1,MXJ,FMAX,FMIN,TR)
C
C Overlay contours in red.
C
      CALL PGSCI(2)
      DO 30 I=1,21
          ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
          IF (MOD(I,5).EQ.0) THEN
              CALL PGSLW(3)
          ELSE
              CALL PGSLW(1)
          END IF
          IF (I.LT.10) THEN
              CALL PGSLS(2)
          ELSE
              CALL PGSLS(1)
          END IF
          CALL PGCONT(F,MXI,MXJ,1,MXI,1,MXJ,ALEV,-1,TR)
   30 CONTINUE
      CALL PGSLS(1)
      CALL PGSLW(1)
C
C Annotate the plot.
C
      CALL PGSCI(2)
      CALL OUTLIN(1,MXI,1,MXJ,TR)
      CALL PGSCI(5)
      CALL PGMTXT('t',1.0,0.0,0.0,'Routines PGGRAY, PGCONT, PGWEDG')
      CALL PGBOX('bcnts',0.0,0,'bcnts',0.0,0)
C
C Draw a wedge.
C
      CALL PGSCH(0.8)
      CALL PGWEDG('BG', 3.0, 4.0, FMAX, FMIN, 'Elevation')
      CALL PGSCH(1.0)
C-----------------------------------------------------------------------
C PGIMAG
C-----------------------------------------------------------------------
C Clear the screen. Set up window and viewport.
C
      CALL PGPAGE
      CALL SETVP
      CALL PGWNAD(-40., 40.,-40., 40.)
      CALL PGSCI(1)
C
C Set up the color map.
C
      CALL PALETT(2)
C
C Draw the map with PGIMAG.  
C
      CALL PGIMAG(F,MXI,MXJ,1,MXI,1,MXJ,FMIN,FMAX,TR)
C
C Overlay contours in white.
C
      CALL PGSCI(1)
      DO 40 I=1,21
          ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
          IF (MOD(I,5).EQ.0) THEN
              CALL PGSLW(3)
          ELSE
              CALL PGSLW(1)
          END IF
          IF (I.LT.10) THEN
              CALL PGSLS(2)
          ELSE
              CALL PGSLS(1)
          END IF
          CALL PGCONT(F,MXI,MXJ,1,MXI,1,MXJ,ALEV,-1,TR)
   40 CONTINUE
      CALL PGSLS(1)
      CALL PGSLW(1)
C
C Annotate the plot.
C
      CALL PGSCI(2)
      CALL OUTLIN(1,MXI,1,MXJ,TR)
      CALL PGSCI(5)
      CALL PGMTXT('t',1.0,0.0,0.0,'Routines PGIMAG, PGCONT, PGWEDG')
      CALL PGBOX('bcnts',0.0,0,'bcnts',0.0,0)
C
C Draw a wedge.
C
      CALL PGSCH(0.8)
      CALL PGWEDG('BI', 3.0, 4.0, FMIN, FMAX, 'Elevation')
      CALL PGSCH(1.0)
C-----------------------------------------------------------------------
C Close the device and exit.
C
      CALL PGEND
      END

      SUBROUTINE OUTLIN(I1,I2,J1,J2,TR)
      INTEGER I1,I2,J1,J2
      REAL TR(6)
C-----------------------------------------------------------------------
C Draw the enclosing rectangle of the subarray to be contoured,
C applying the transformation TR.
C
C For a contour map, the corners are (I1,J1) and (I2,J2); for
C a gray-scale map, they are (I1-0.5,J1-0.5), (I2+0.5, J2+0.5).
C-----------------------------------------------------------------------
      INTEGER K
      REAL XW(5), YW(5), T
C
      XW(1) = I1
      YW(1) = J1
      XW(2) = I1
      YW(2) = J2
      XW(3) = I2
      YW(3) = J2
      XW(4) = I2
      YW(4) = J1
      XW(5) = I1
      YW(5) = J1
      DO 10 K=1,5
          T = XW(K)
          XW(K) = TR(1) + TR(2)*T + TR(3)*YW(K)
          YW(K) = TR(4) + TR(5)*T + TR(6)*YW(K)
   10 CONTINUE
      CALL PGLINE(5,XW,YW)
      END

      SUBROUTINE PALETT(TYPE)
C-----------------------------------------------------------------------
C Set a "palette" of colors in the range of color indices used by
C PGIMAG.
C-----------------------------------------------------------------------
      INTEGER TYPE
C
      REAL GL(2), GR(2), GG(2), GB(2)
      REAL RL(9), RR(9), RG(9), RB(9)
      REAL HL(5), HR(5), HG(5), HB(5)
      REAL WL(10), WR(10), WG(10), WB(10)
      REAL AL(20), AR(20), AG(20), AB(20)
      REAL TL(4), TR(4), TG(4), TB(4)
C
      DATA GL /0.0, 1.0/
      DATA GR /0.0, 1.0/
      DATA GG /0.0, 1.0/
      DATA GB /0.0, 1.0/
C
      DATA RL /-0.5, 0.0, 0.17, 0.33, 0.50, 0.67, 0.83, 1.0, 1.7/
      DATA RR / 0.0, 0.0,  0.0,  0.0,  0.6,  1.0,  1.0, 1.0, 1.0/
      DATA RG / 0.0, 0.0,  0.0,  1.0,  1.0,  1.0,  0.6, 0.0, 1.0/
      DATA RB / 0.0, 0.3,  0.8,  1.0,  0.3,  0.0,  0.0, 0.0, 1.0/
C
      DATA HL /0.0, 0.2, 0.4, 0.6, 1.0/
      DATA HR /0.0, 0.5, 1.0, 1.0, 1.0/
      DATA HG /0.0, 0.0, 0.5, 1.0, 1.0/
      DATA HB /0.0, 0.0, 0.0, 0.3, 1.0/
C
      DATA WL /0.0, 0.5, 0.5, 0.7, 0.7, 0.85, 0.85, 0.95, 0.95, 1.0/
      DATA WR /0.0, 1.0, 0.0, 0.0, 0.3,  0.8,  0.3,  1.0,  1.0, 1.0/
      DATA WG /0.0, 0.5, 0.4, 1.0, 0.0,  0.0,  0.2,  0.7,  1.0, 1.0/
      DATA WB /0.0, 0.0, 0.0, 0.0, 0.4,  1.0,  0.0,  0.0, 0.95, 1.0/
C
      DATA AL /0.0, 0.1, 0.1, 0.2, 0.2, 0.3, 0.3, 0.4, 0.4, 0.5,
     :         0.5, 0.6, 0.6, 0.7, 0.7, 0.8, 0.8, 0.9, 0.9, 1.0/
      DATA AR /0.0, 0.0, 0.3, 0.3, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0,
     :         0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/
      DATA AG /0.0, 0.0, 0.3, 0.3, 0.0, 0.0, 0.0, 0.0, 0.8, 0.8,
     :         0.6, 0.6, 1.0, 1.0, 1.0, 1.0, 0.8, 0.8, 0.0, 0.0/
      DATA AB /0.0, 0.0, 0.3, 0.3, 0.7, 0.7, 0.7, 0.7, 0.9, 0.9,
     :         0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/
C
      DATA TL /0.0, 0.5, 0.5, 1.0/
      DATA TR /0.2, 0.6, 0.6, 1.0/
      DATA TG /0.0, 0.0, 0.5, 1.0/
      DATA TB /1.0, 0.0, 0.0, 0.0/
C
      IF (TYPE.EQ.1) THEN
C        -- gray scale
         CALL PGCTAB(GL, GR, GG, GB, 2, 1.0, 0.5)
      ELSE IF (TYPE.EQ.2) THEN
C        -- rainbow
         CALL PGCTAB(RL, RR, RG, RB, 9, 1.0, 0.5)
      ELSE IF (TYPE.EQ.3) THEN
C        -- heat
         CALL PGCTAB(HL, HR, HG, HB, 5, 1.0, 0.5)
      ELSE IF (TYPE.EQ.4) THEN
C        -- weird IRAF
         CALL PGCTAB(WL, WR, WG, WB, 10, 1.0, 0.5)
      ELSE IF (TYPE.EQ.5) THEN
C        -- AIPS
         CALL PGCTAB(AL, AR, AG, AB, 20, 1.0, 0.5)
      ELSE IF (TYPE.EQ.6) THEN
C        -- TJP
         CALL PGCTAB(TL, TR, TG, TB, 4, 1.0, 0.5)
      END IF
      END

      SUBROUTINE SETVP
C-----------------------------------------------------------------------
C Set the viewport, allowing margins around the edge for annotation.
C (This is similar in effect to PGVSTD, but has different margins.)
C The routine determines the view-surface size and allocates margins
C as fractions of the minimum of width and height.
C-----------------------------------------------------------------------
      REAL D, VPX1, VPX2, VPY1, VPY2
C
      CALL PGSVP(0.0, 1.0, 0.0, 1.0)
      CALL PGQVP(1, VPX1, VPX2, VPY1, VPY2)
      D = MIN(VPX2-VPX1, VPY2-VPY1)/40.0
      VPX1 = VPX1 + 5.0*D
      VPX2 = VPX2 - 2.0*D
      VPY1 = VPY1 + 8.0*D
      VPY2 = VPY2 - 2.0*D
      CALL PGVSIZ(VPX1, VPX2, VPY1, VPY2)
      END
